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Abstract. Polyhomeostatic adaption occurs when evolving systems try to achieve a 
target distribution function for certain dynamical parameters, a generalization of the 
notion of homeostasis. Here we consider a single rate encoding leaky integrator neuron 
model driven by white noise, adapting slowly its internal parameters, the threshold 
and the gain, in order to achieve a given target distribution for its time-average firing 
rate. For the case of sparse encoding, when the target firing-rated distribution is 
bimodal, we observe the occurrence of spontaneous quasi-periodic adaptive oscillations 
resulting from fast transition between two quasi-stationary attractors. We interpret 
this behavior as self-organized stochastic tipping, with noise driving the escape from 
the quasi-stationary attractors. 
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1. Introduction 

Self-regulation plays an important role in biological and technical systems. 
Homeostatically regulated steady states are a precondition to life, examples being the 
concentration of blood glucose controlled by insulin [1] and glucagon, the pH value 
of blood [21 E] and the body temperature [2], which are all autoregulated in order to 
maintain stable conditions. Further examples are the concentration of ions, proteins and 
transmitters in the brain, their respective levels are all self regulated [5]. Furthermore, 
homeostasis is implemented and can be found in technical systems, for example in 
microrobotic swarms [5J. Adaption typically introduces a slow time scale into the 
dynamical system [7], a process also denoted meta learning, a central notion in the 
context of neuromodulation [H] and emotional control [5]. The resulting dynamical 
system then has both fast and slow variables and critical transitions in the form of 
tipping processes may occur [TO] . 

Classical homeostasis involves the regulation of a scalar quantity, like the body 
temperature. More complex forms of homeostasis are however also important in the 
realm of life. For example, an animal may want to achieve a certain time averaged 
distribution of behaviors, like foraging, resting and engaging socially, over the course 
of several days. This kind of adaptive behavior has been termed polyhomeostasis 
[TT1 [12] . It occurs when a dynamical system tries to achieve, via the continuous 
adaption of slow variables, a given target distribution for the time-averaged activity 
of a subset of fast variables. Polyhomeostatically adapting systems are typically slow- 
fast dynamical systems and their dynamical behavior can tip spontaneously from one 
state into another. For a network of rate-encoding neurons tipping transitions from 
laminar flow to intermittent chaotic bursts of activities have been observed [TT1 IT2] . 

Tipping transitions can occur both in adaptive and in driven systems. Potential 
tipping scenarios are currently discussed intensively in the context of climate research 
[131 [H], they may be related to a slow driving of external parameters [15J, to noisy input 
inducing a stochastic escape from a local attractor [IS1EI7J, or through a dynamical effect 
when the rate of change of a control parameter reaches a certain threshold [13J. 

Here we study the phenomenon of self-organized tipping for a polyhomeostatic 
adapting system driven by a steady-state stochastic input. We examine a previously 
proposed model [HI [19] for regulating the firing rate distribution of individual neurons 
based on information-theoretical principles. This type of model has been studied 
previously for the case of discrete time systems and unimodal target firing rate 
distributions [TTl [T2"] . Here we examine the case of continuous time and bimodal 
target distribution functions, corresponding to sparse coding. For bimodal firing rate 
distributions the neural activity tends to switch in between states close to minimal 
and maximal activity. Similar bimodal activity states are observed also in many other 
systems, e.g. dynamical gene regulation networks [20]. We find that bimodal target 
distributions may lead to self-organized bistability within a certain range of parameters. 

We consider a single leaky integrator neuron with noisy input and a sigmoidal 
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Figure 1. The transfer function g(x), see Eq. ([2]), for thresholds b — 2 (red lines) and 
= 3 (green lines) and various gains a: 1 (dotted), 3 (dashed), 9 (solid). 

transfer function having two degrees of freedom. To achieve a special behavior - here 
the temporal output distribution of the firing rate - we use polyhomeostasis to change 
the intrinsic parameters which are directly influencing the mapping of the membrane 
potential to the firing rate in order to obtain a specific output distribution. We derive 
these parameter changing rules using stochastic adaption and show that two degrees of 
freedom already result in a good behavior approximation, for most of the parameters 
studied. For bimodal adaption target distributions we observe self-organized and quasi 
periodic stochastic tipping in between two quasi-stationary attractors resulting from 
competing adaption gradients. 

2. Model 

Biological neurons integrate incoming signals and emit an axon potential, a spike, 
whenever the membrane potential has reached a certain threshold. The membrane 
potential then returns, after a short refractory period, rapidly to its resting value. This 
behavior can be captured using spiking integrate and fire neural models [21]. In many 
circumstances the firing rate, the number of spikes per unit time, is important and rate 
encoding neural models can be used [22]. Here we consider a single rate-encoding leaky 
integrator driven by white noise £(t), 

x(t) = -Tx(t) + £(*), (£(*)£(0> = QS(t - , (1) 

where x > is the membrane potential and T > the relaxation rate. The firing 
rate y(t) G [0, 1] is a nonlinear function of the membrane potential x(t), which we have 
selected as 

y(t) = g(x(t)), 0( X ) = __JL^, (2) 

where a > is the gain and b is the threshold. The polynomial transfer function ^ 
maps the membrane potential x G [—00,00] to the normalized firing rate y G [0,1] 
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which approaches zero and unity for small and large membrane potentials respectively, 
compare Fig. [TJ The slope of g(x) is a/4 at the threshold b. 

Usually the intrinsic parameters of the transfer function ^2§, a and b are taken as 
given by some a priori considerations. Here we will consider them to be slow variables, 
a = a(t) and b = b(t), adapting slowly such that a target dynamical behavior is 
approached on the average for the firing rate y{t). The stochastic driving G [Hi, 5 2 ] 
in (jl]) is simulated through white noise plateaus: The values are generated according to a 
uniform probability distribution (white noise), but they remain constant for short time 
intervals on the order of unity. The membrane potential averages the input driving 
noise, due to the leak rate r in (jl]), its distribution function p(x) having a mean 
p p « (Hi + S 2 )/(2r) and variance a 2 p « (S 2 - Si)/(2r). 



2.1. Polyhomeo static Adaption 
The firing-rate statistics is given by 

p(z) = - 6(z-y(t))dt, p(z)dz = l, (3) 

1 J to JO 

where the length T of the sliding observation window is substantially larger than 
the relaxation rate 1/T. The firing-rate distribution p(z) is an important quantity 
characterizing the information processing capability of biological and artificial neurons. 
No information is encoded for a constant firing rate, the next value is always exactly 
the same as before, so no new information is transferred. One may assume that a 
certain distribution q(y) of firing rates may constitute an optimal working regime. 
Possible functional dependencies for q(y) can be derived by information-theoretical 
considerations, e.g. maximizing information entropy, as discussed further below. 

Considering a given target firing-rate distribution q(y), the closeness of the actual 
firing-rate distribution p(y) is measured by the Kullback-Leibler divergence (KL- 
divergenceO D KL , their relative entropy, [7]: 

Dkl(p, q)= I dyp(y) In ?f\, D KL (p, q) > . (4) 

The Kullback-Leibler divergence is positive definite and vanishes only when the two 
distribution coincide. The KL-divergence is generically not symmetric but becomes 
symmetric in the limiting case of similar distributions p and q, becoming equivalent 
in this limit to the \ 2 test [7]. Our aim is now to rewrite Q as an integral over the 
membrane potential x, using 

P\X) 

P(y)dy = p(x)dx, p(y) = , (5) 

dy/dx 

where p(x) is the membrane potential distribution. Using y = g(x) and Eqs. Q and 
(j5|, we obtain 



ODkl 
06 



dx p(x) 



1 dg' _ q[_dg_ 
V 09 q 09 



= l dx P( x )^jj ( 6 ) 
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Figure 2. Target distribution q(y), see ([9]), with some selected parameters Ai and A2. 
The target firing-rated distributions is bimodal for A2 > 0. 



for the derivative of the Kullback-Leibler divergence with respect to the intrinsic 
parameters 9 = a, b of the transfer function g(x), see (]2j. 

We consider now the case that the system does not dispose of prior information 
about the distribution of input stimuli and the thereby resulting distribution of 
membrane potential p(x). The best strategy to minimize the Kullback-Leibler is then to 
minimize the individual terms of the integral (|6| through the stochastic adaption rules 

mm 

Tt = - ee w 9 = a ' b (7) 

for the intrinsic parameters of the transfer function g(x), where the eg are appropriate 
small adaption rates. 

2.2. Target Firing-Rate Distribution 

In order to evaluate Q, respectively Eq. ([6]), we need to specify the target firing rate 
distribution q(y). For this purpose we use information-theoretical considerations. 

Given a continuous probability distribution function q its Shannon entropy H(q) 
can be defined as 

H{q) = - [ dyq(y)\nq(y) (8) 



Among all the real-valued distributions with specified mean \x and standard deviation 
a the Gaussian distribution [7] 

q(y) oc exp I 1 oc exp (\ t y + \ 2 y 2 ) (9) 
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Figure 3. Typical time series for a mono-modal target distribution q(y) with 
Ai = —10, A2 = 0, compare Fig. [2j Plotted are the membrane potential x (solid 
blue line, upper panel), the threshold b (dashed red line, upper panel), the gain a 
(solid green line, middle panel) and the firing rate y (solid black line, lower panel). 
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has maximal information entropy, with \i 
obtained using variational calculus, 



— Ai/ (2A2) and 2cr 2 = — I/A2, which is easily 



= 5 



H{q) + X 1 J dyyq{y) + X 2 



dyy 2 Q(y) 



where (— Ai) and (— A2) are the respective Lagrange parameters. In Fig. [2] examples for 
q(y) are illustrated for several values of Ai and A 2 . The support of the target firing rates 
is compact, y G [0, 1], and both negative and positive Ai and A 2 can be considered. The 
normalization factor dy q(y) cancels out in (6J), since only ratios are involved. 

For positive A 2 > and Xi w — A 2 one obtains bimodal target distributions. This 
is an interesting case, since sparse coding, which is realized when only a minority of 
neurons of a given network is active, and a majority is inactive [23], is characterized by 
a skewed bimodal distribution. 



2.3. Stochastic Adaption Rules 
From ([9]) and ^ we find the relations 



q(y) 



Ai + 2X 2 y, 



dx 



ag (1 - 9) 
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and 

= (x - b)g (1 - g) , ^ = -05 (1 - o) , 
which we can use to evaluate the stochastic adaption rules ^ as 



and 



da 
dt 

db 
dt 



+ (x-6) l-2y + {\ 1 + 2\ 2 y){l-y)y 



-a(l-2y + (X 1 + 2X 2 y) (l-y)y 



(10) 



These two adaption rules will lead to an adaption of the time-averaged firing rate 
distribution p(y) towards the target distribution q{x) whenever the adaption time-scales 
1/eg are substantially larger than the time constants of the neural dynamics, which in 
turn are determined by the time scale of the incoming stimuli and by the leak-rate T in 
Q. ^ 

The transfer function g(x) contains only two free parameters, the gain a and the 
threshold b. Perfect adaption p(y) = q(y), for all y G [0, 1] can hence not be expected. 
The system tries to minimize the Kullback-Leibler divergence by adapting the available 
degrees of freedom, which are just two in our case. 



2.4- Numerical method 

The equations Q, ( [To] ) and (11) form a set of first order differential equations with 
respect to time. We solve them numerically using the Euler method with one evaluation 
per time step. The random white noise is generated through a pseudo-random number 
generator with a uniform distribution. The values for the leak T, the time step At and 
the learning rates e a and % are shown in the corresponding figures. 



3. Results 



We performed a series of simulations with the aim to study two issues. Polyhomeostatic 
adaption had been studied previously for the case of discrete time systems [TH [11] , here 
we examine the case of continuous time. The case of a bimodal target distribution is, 
in addition, highly interesting, as it confronts the system with a dilemma. The transfer 
function g(x), compare Fig. [TJ is strictly monotonic. The distribution of the membrane 
potential p(x) is hence mono-modal. There is no easy way for the adapting neuron to 
achieve, as a steady state time-average, a bimodal output firing rate distribution p(y). 
The question then is whether the system will find a way out of this dilemma through 
spontaneous behavioral changes. 



3.1. Target Distribution Approximation 

For most simulations we used, if not stated otherwise, T = 1 for the leak rate and 
At = 10 _1 for the integration time step. A typical time series is given in Fig. |3j Note 
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Figure 4. Target distribution q (bars) vs. achieved distribution p (points) for different 
distributions. Ai and A2 are given in each diagram. At = 10 , t niax = 10 8 , 
e a = e 6 = 10- 2 , T = 1, 5 =[0,10]. 



that the adaption of the intrinsic parameters a and b takes place on a slower time 
scale than the one of the primary dynamic variables, x and y, as typical for a slow-fast 
dynamical system. 

Applying moderate to small learning rates e a = et, < 0.01 the neuron's firing rate y 
approximates various types target distributions q quite well. In Fig. [4] the achieved and 
the respective target firing rated distributions are compared. The respective relative 
entropies are well minimized and presented in Table [T] Strictly speaking the stochastic 
adaption rules (10) and (11) are equivalent to approximating the firing-rate statistics 
Q, which is a time-averaged quantity, towards the target distribution function q(y) 
only in the limit of very small adaption rates, e a and e&. Small but finite values for 
the adaption rates, as used in our simulations, correspond to to a trailing averaging 
procedure over a limited time interval, and the value of Kullback-Leibler divergence 
achieved hence depend weakly on the actual values used for the learning rates. 

For very high learning rates, e& 3> 0.1, the threshold b follows the membrane 
potential x nearly instantaneously, both variables become highly correlated. Therefore 
the firing rate distribution p cannot approximate the target distribution q any more, in 
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fact the resulting Kullback-Leibler divergence is then very high. The tipping in dynamic 
behavior as a function of adaption rate amplitude is typical for a rate-induced tipping 
transition [13] . 

3.2. Gain- Threshold Phase Diagram 

Due to the sigmoidal shape of the transfer function, several target distributions lead 
to specific fingerprints in the gain-threshold phase diagram which we present in Fig. [5j 
The threshold, for example, for a left (right) dominant target distribution is high (low) 
and is therefore sensitive to the mean /i = — Ai/(2A2) of q{y). Small gains a result in 
quite flat transfer functions g(x), compare Fig. [TJ mapping the membrane potentials to 
similar firing rates y. High gains a discriminate, relative to the threshold b, on the other 
side between high and low membrane potentials. The gain is therefore smaller for hill 
shaped and flat target distributions, as compared to the left and right dominant target 
distributions (e.g. Ai = —20, A2 = +20) for which intermediate values are suppressed. 

Left (right) dominant target distributions (compare Table [T]) correspond directly to 
high (low) transfer function thresholds. Uniform, hill and other not unilateral dominant 
target distributions lead to intermediate transfer function thresholds with a wide variety 
of the transfer function gains. For symmetrical target distributions from hill shaped to 
diametrical shaped there is a transition from low to high gains. 

3.3. Self- Organized Stochastic Escape 

While the left or right dominant target distributions are easily approximated due to 
the sigmoidal shape of the transfer function, the bimodal left and right dominant target 
distributions puts the system in dilemma: Since intermediate values are to be suppressed 
the transfer function gain a cannot be too small. Because of this there exists at least two 
quasi-stationary fixed points, one for the left, one for the right part of the distribution. 

For zero or small learnings rates e a = e& « the system is trapped in a single local 
fixed point. Only the left or right part of the target distribution is then approximated, 

Table 1. The relative entropies Dkl Q of various target distributions (see Fig. [2]) 
compared to the corresponding achieved distribution, compare Fig. [4j 
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Figure 5. Phase diagram. Plotted are the gain a(t) and the threshold b(t) of the 
transfer function for various target distributions (Ai and A2 given in the legend). The 
respective target and achieved firing rate distributions are given in Fig. [4] 



the Kullback-Leibler divergence is not well minimized. 

Increasing the learning rates e a = e& allows the system to escape stochastically from 
the respective local fixed points: The transfer function threshold b conquers the local 
gradient and moves to the other fixed point and back (compare Fig. [6]). In the long-term 
observation the system therefore approximates the left and the right part of the target 
distribution and hence minimizes the relative entropy, compare Table [2j These tipping 
transitions between the two quasi-stationary fixed points are illustrated in Fig. [7j which 
shows a typical time series for a skewed target distribution. Note that there are two 
fixed points for the gain and threshold and a direct correspondence to the periods of 
high and low firing rates y(t). 



Table 2. Relative entropies Dkl @ for the left-skewed target distribution 
(Ai = —20, A2 = 18.5) relative to the achieved distribution for various learning rates 
e a and ej,, compare Fig [6] 
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Dkl 



0.306 0.295 0.293 0.289 



0.283 0.154 



0.109 




Figure 6. Stochastic escape: Phase diagram of the transfer function gain vs. transfer 
function threshold for a convex left skewed target distribution with various learning 
rates (e a = ej given in the legend). At = 10 _1 , V = 1, Ai = —20, A 2 = 18.5. 



Very low learning rates e a , e?,, lead to deep and big basins of attraction for the 
respective fixed points, while on the other hand high learning rates result in the closely 
following of the threshold to the membrane potential which prohibits reaching the target 
distribution. This mechanism is reminiscent to the case of Langevin dynamics in a 
double- well potential [23], where a stochastically driven particle may switch forth and 
back between two local minima [7]. The switching time is controlled for the double- well 
problem by the Kramer's escape rate, which depends exponentially on the potential 
barrier height. It is difficult to formulate a quantitative mapping to the double- well 
problem, the local attractors visible in Figs. [6] and [7], and the effective barriers in between 
them, are self-organized structures. Note that the strength Q of the noise term is 
constant and influences the transition rate only weakly, due to the continuous adaption 
of the transfer function, via (10) and (11), to the average strength of the stochastic 
driving. 



4. Discussion 



We showed that polyhomeostatic adaption of continuous-time leaky integrator leads to 
desired firing rate distributions. We also run further simulations using white noise and 
Gaussian noise input and replace the transfer function by other qualitatively different 
(but still sigmoidal) functions, see Appendix. It turns out that the polyhomeostatic 
adaption as well as the self-organized stochastic escape are quite robust principles. 
However, the quality of the approximation (as seen by visual overlapping) and the 
value of the Kullback-Leibler divergence depend on the learnings rates, but also on the 
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Figure 7. Time series: membrane potential x, transfer function threshold b (clashed), 
transfer function gain a and firing rate y. At = 10 _1 , e a = et, = 10 -1 , r = 1, Ai = —20, 
A 2 = 18.5. 

input distribution and the input's strength. 

The stochastic tipping as a function of adaption rates has a close relation to the 
phenomenon of stochastic escape. The strength of the driving input noise is constant, 
but its influence is averaged out for very low adaption rates. Stochastic escape from 
one local attractor to another is not possible. The stochasticity of the input becomes 
relevant for intermediate values of adaption rates and stochastic transitions between 
the two quasistationary attractors are most frequent. Finally, for very large adaption 
rates, the system tips into another dynamical state, tracking the stochastic input signal 
nearly instantaneously. This sequence of behaviors is self organized and reached from 
any initial state. 
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Figure AI. The transfer function g(x), see Eq. (A.l), for thresholds 6 = 2 (red lines) 



and b = 3 (green lines) and various gains a: 1/3 (dotted), 3 (dashed), 9 (solid). No 
inflection point is present for exponents ab < 1. 



Appendix A. Polynomial transfer function 



The polyhomeostatic adaption of the system is not changing qualitatively by replacing 
the transfer function g. Instead it turns out that the system is robust against changing 
the transfer function as long as it remains sigmoidal. We also applied a transfer function 

9{x) = , {X 'T. ~ , (A.l) 



(x/b) ab + 1 ' 

with a polynomial decay to g(0) = 0, which limits the membrane potential x > to 
be non-negative. It turns out that the shape of the target distribution q is also well 
approximated using this transfer function. Also stochastic escape from one fixed point 
to another and back can be observed as for some target distributions two fixed points 
are necessary. 

The transfer function has an inflection point for exponents ab > 0; it is absent for 



ab < 1, compare Fig. |A1| The transfer function g behaves as 

(x/bf b 

2 H - ~^^ J ~ 



x <^b 
b) x ~ b 
x^>b 



g(x) w { \ + \a{x-b) x^b . (A. 2) 

1 - (b/x) 

The slope is a/4 which approaches zero and unity for small and large membrane 
potentials respectively. 



From Eq. (A.l) we find the relations 

ab 
9)9- 



dg_ 
dx 
dg_ 
da 



1 



,r 



1 - g) gb\\i 



x 



dg_ 
db 



;i -g) ga 



In 



which we can use to evaluate the stochastic adaption rules ([7]) as 



da 
di 



b\n{x/b) l-2y+{\ 1 + 2\ 2 y) (l-y)y 



(A.3) 
(A.4) 

(A.5) 
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and 



db 
dt 



\n(x/b) — 1 



2y + {\ 1 + 2X 2 y){l-y)y 



(A.6) 



Applying this transfer function g it turns out that the target distribution is well 
approximated also in this case, even though the membrane potential is restricted to 
non-negative numbers. Tab. Al lists the well minimized Kullback-Leibler divergences 
for several target distributions. 

We conclude that the stochastic adaption rules are therefore generic and 
qualitatively independent on the concrete realization of the transfer function. However, 
quantitatively the resulting relative entropies depend on the choice of the transfer 
function which also influences the optimal adaption rates e a and e&. 



Appendix A.l. Self-Organized Stochastic Escape 

For the non-symmetric convex target distribution (Ai = —20, A2 = 19) there are two 
fixed points. Since the target distribution cannot be well approximated by only one 
fixed point the system escapes stochastically from one to the other and back with a 
certain period, compare Fig. A2 For small learning rates e a = e& ;$ 0.01 the system in 



trapped in only one fixed point. The relative entropy therefore is not well minimized. 
For intermediate learning rates 0.01 ^ e a = €b ~ 0.04 the perturbation is high 



enough to stochastically escape from that fixed point and approach another one. Fig. |A2 



Table Al. Relative entropies of various target distributions compared to the 
corresponding achieved distribution (e a = = 10~ 2 , bins = 100). 



Ai 


A 2 


shape 


Dkl 








uniform 


0.060131 


-10 





left dominant 


0.069351 


+10 





right dominant 


0.114578 


-10 


+10 


left/right dominant 


0.051811 


+20 


-20 


hill 


0.148098 


-20 


+20 


left/right, symmetric 


0.189217 


-20 


+19 


left /right, left skewed 


0.063934 


-20 


+18.5 


left /right, left skewed 


0.261215 



Table A2. Relative entropies of the left-skewed target distribution (Ai = —20, 
A2 = 19) compared to the achieved distribution for various learning rates e a and 
£fc. Note that the Kullback-Leibler divergence is not minimized for Cb ~ 0.05 due to 
the fast correlation of the membrane potential and the transfer function threshold. 



10- 



10- 



0.01 0.03 0.04 0.05 



0.06 



D 



KL 



0.376 0.368 0.064 0.043 0.017 1.892 1.591 
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gain a 



Figure A2. Left: Stochastic escape: Phase diagram of the transfer function gain vs. 
transfer function threshold for a convex left skewed target distribution with various 
learning rates (e a and eb given in the legend). At = 10~ 2 , T = 0.1, Ai = —20, A 2 = 19. 
Right: Time series: membrane potential x, transfer function threshold b (dashed), 
transfer function gain a and firing rate y. At = 10 _1 , e a = €b = 10~ 2 , T = 0.1, 
Ai = -20, A 2 = 19. 



shows a typical time series for this tipping. This has also an effect on the relative entropy 
which therefore is even smaller than without stochastic escape (see Tab. A2). 

For high learning rates e a = e& ^ 0.05 the system's behavior changes: the transfer 
function is close to a Heaviside step function and the threshold follows the membrane 
potential quickly. In that state the achieved distribution is not close to the target 
distribution, therefore the relative entropy is not minimized anymore (see Tab. A2). 



